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A  wide  variety  of  line-drawing  algorithms  can  be  derived  by  applying  program  transformations  to  a 
simple,  obviously  correct  algorithm.  The  transformations  increase  the  algorithm's  performance  and 
eliminate  the  need  for  floating-point  computations.  Two  familiar  algorithms  are  derived  in  this 
way:  Urcsenham's  algorithm  and  the  digital  differential  analyzer  (DDA).  The  transformations  are 
then  used  to  derive  several  highly  parallel  variants  of  Bresenham's  algorithm,  designed  for  use  on 
displays  that  can  generate  more  than  one  pixel  at  a  time.  ITic  treatment  shows  a  complete, 
extended  example  of  the  practical  use  of  program  transformations.  Moreover,  the  transformations 
derive  Bresenham’s  algorithm  without  recourse  to  complex  geometric  arguments. 
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1.  Introduction 

Many  computer  graphics  devices  use  "line-drawing  algorithms"  to  approximate  the  appearance 
of  straight  lines  on  devices  that  can  only  produce  dots  on  a  discrete  grid.  Incremental  pen  plotters 
that  move  a  pen  in  small  steps  arc  common  devices  that  require  such  a  line-generation  algorithm. 
Point-plotting  CR  T  displays  use  the  algorithms  to  approximate  straight  lines  on  interactive  graphics 
displays.  More  recently,  frame-buffer  raster-scan  displays  use  these  algorithms  to  identify  the 
picture  elements  (pixels)  that  should  be  illuminated  to  display  a  line. 

Simplicity  and  speed  the  the  key  design  criteria  for  line-drawing  algorithms,  because  the 
computations  arc  often  implemented  in  hardware  in  order  to  achieve  high  line-generation  speeds.  It 
appears  that  the  early  popularity  of  the  binary  rate  multiplier  (URM)  was  due  entirely  to  simplicity, 
for  it  generates  rather  poor  approximations  to  straight  lines  Ihc  digital  differential  analyzer  (DDA) 
generates  better  approximations  to  die  true  line,  but  requires  an  iterative  loop  that  may  average 
almost  two  cycles  to  generate  each  point.  An  algorithm  devised  by  J.H.  Hresenham  [l|  dominates 
the  DDA:  it  generates  the  optimal  line,  in  a  sense  of  optimal  described  below;  it  requires  only 
integer  additions  and  subtractions;  and  it  generates  one  output  point  for  each  iteration  of  the  inner 
loop. 

To  achieve  very  high  line-generation  speeds,  we  need  algorithms  that  can  determine  the 
location  of  several  points  on  a  line  in  parallel.  None  of  the  current  line-drawing  techniques  is 
suitable,  as  they  trace  out  the  line  sequentially,  one  point  at  a  time.  Parallel  algorithms  have  several 
applications,  chiefly  in  raster-scanned  systems  that  can  write  more  than  one  pixel  at  a  time  into  the 
image.  Ihc  "8x8  frame  buffer  display"  [2],  which  can  in  one  memory  cycle  write  a  square  region  8 
pixels  on  a  side  located  anywhere  on  the  screen,  motivated  the  investigation  of  parallel  algorithms. 

Ibis  paper  shows  how  simple  program  transformations  are  used  to  derive  all  of  these 
algorithms,  starting  from  obviously  correct  algorithms  based  on  simple  analytic  geometry.  These 
transformations  assure  us  th.it  the  more  efficient  but  more  complex  algorithms  are  correct,  because 
they  have  been  derived  by  correct  transformations  from  a  correct  algorithm. 


2.  Line-drawing  preliminaries 

'The  line-drawing  problem  is  to  dcterr.ine  a  set  of  pixel  coordinates  (v,  y).  where  v  and  y  arc 
integers,  that  closely  approximates  the  line  from  die  point  (0,0)  to  the  point  (d\.  dy).  for  integer 

values  of  dx  and  dv.  The  assumption  that  one  line  endpoint  is  at  (he  origin  loses  no  generality, 

because  lines  with  other  origins  arc  simply  translations  of  the  line  with  origin  (0.0).  Additionally, 

lines  arc  restricted  to  the  first  octant:  0  <  dy  <  dx.  Again,  this  assumption  lose-,  no  generality, 

because  an  arbitrary  lino  can  be  generated  by  transposing  the  canonical  line  or  by  reflecting  it  about 
one  of  the  principal  axes. 
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Ihc  objective  of  a  line-drawing  algorithm  is  to  enumerate  those  pixels  that  lie  closest  to  the 
true  line,  the  mathematical  line  from  (0.0)  to  (</x  dy).  l-'igurc  1  illustrates  a  typical  line,  showing 
with  circles  the  pixels  that  correspond  to  spots  illuminated  by  a  CRT  beam  on  a  raster  display  or  to 
the  swath  of  a  plotter  pen.  Notice  that  integral  values  of  coordinates  locate  pixel  centers. 


Figure  I.  The  line  from  (0.  0)  to  (8.  5).  Small  dots  represent  pixel 
centers.  The  solid  line  represents  the  "true"  line.  Circles  show  the 
pixels  that  are  illuminated  to  display  the  optima!  line. 

The  optimal  line  will  illuminate  exactly  one  pixel  in  each  vertical  column.  Ihis  assumption 
minimizes  variations  in  pixel  spacing  that  make  lines  appear  to  vary  in  width  or  brightness,  lire 
assumption  depends  on  the  fact  that  the  line's  extent  in  x  exceeds  its  extent  in  y. 

The  line-drawing  algorithm  must  compute,  for  each  integer  xf  the  coordinate  y  of  the  pixel 
that  should  be  illuminated.  The  coordinate  y  t  of  the  true  line  is  simply  yt  =  (dy/ilx)Xj. 
Illuminating  a  pixel  centered  at  y{  introduces  an  error  ev  =  j(- =  y[-idy/Jx)xr  measured  along 
the  y  axis.  The  error  c ^  measured  perpendicular  to  the  line  can  be  determined  using  similar 
triangles  (Figure  2):  =  (dx/V (dx2  y  dy))ev  Thus,  for  any  given  line,  e is  simply  a  constant 

times  er  Consequently,  determining  y(  by  minimizing  the  error  r,,  will  identify  the  pixel  that  is 
closest  to  the  line,  using  cither  vertical  or  perpendicular  distance  measures. 

The  errors  can  be  minimized  if  y(  is  computed  by  rounding  yf:  y-  =  mun iKy,).  or  equivalently. 
y  -  trunciyl+ 1/2)  =  Ly;  +  1/2J.  (Recall  that  the  floor  function.  LxJ.  denotes  die  greatest  integer 
less  than  or  equal  to  x.)  With  this  choice,  cy  =  Ly^-t- l/2J-v.,  so  -  1/2  <  cy  <  1/2.  Lines  with 
this  error  behavior  arc  said  to  be  optimal,  in  the  sense  that  each  pixel  illuminated  is  within  one-half 
unit  of  the  true  line.  Optimality  thus  requires  that  a  single  pixel  be  illuminated  in  each  column  and 
dial  the  pixel  be  the  one  closest  to  the  true  line 
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dy 


dx 

Figure  2.  Illustration  of  the  relationship  between  the  vertical 
distance  ey  and  the  perpendicular  distance  ey 


3.  Derivation  of  the  Rresenham  algorithm 

The  minimum-error  formulation  of  the  optimal  line  leads  directly  to  a  simple  algorithm  that 
enumerates  all  the  points  on  the  optimal  line,  which  can  be  expressed  in  a  PasCAI -like  language: 

Al:  var  yt:  exactreal;  dx,  dy,  xi,  yi:  integer; 

for  xi  :  =  0  to  dx  do  begin 
yt  :  =  [dy/dx|*xi; 
yi  :  =  trunc(yt  +  [l/2]); 
display(xi.yi) 
end 

Although  this  procedure  is  expressed  using  programming-language  constructs,  it  requires  that  precise 
real  arithmetic  is  used;  "floating-point"  approximations  are  not  permitted.  To  emphasize  this 
precise  arithmetic,  variables  that  use  it  arc  declared  to  have  type  exactreal.  Square  brackets  enclose 
expressions  whose  values  do  not  change  during  iteration  of  the  loop;  these  expressions  can  be 
computed  only  once,  before  the  loop  is  entered,  and  saved  in  temporary  variables.  We  shall  also 
maintain  that  multiplications  by  a  power  of  two  do  not  require  multiplication  operations,  but  can  be 
achieved  by  addition  or  arithmetic  shifting. 

1.  Incremental  transformation.  Ihc  next  version  of  the  algorithm  is  derived  from  Al  by 
observing  that  yt  can  be  calculated  incrementally  by  adding  the  quantity  (dy/dx)  on  each  iteration. 
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A2:  var  yt:  exactreal:  dx.  dy,  xi,  yi:  integer; 

yt  :=  0; 

for  xi  :  =  0  to  dx  do  begin 

yi  :=  trunc(yt  +  [l/2]);  (•  assen  yt  =  (dy/dx)xi  •) 

display(xi.yi); 

yt  :=  yt+[dy/dx] 

end 

2.  Substitution  of  variable  (simple).  A  simple  transformation  substitutes 

ys  =  >,+  1/2  (l) 

A3:  var  ys:  exactreal:  dx,  dy,  xi,  yi:  integer; 

ys  :=  1/2; 

for  xi  :=  0  to  dx  do  begin 

yi  :=  trunc(ys);  (*  assen  ys  =  (dy/dx)xi  +  l/2  =  yt  +  1/2  •) 

display(xi.yi); 

ys  :=  ys  +  [dy/dx] 

end 

3.  Substitution  of  variable  (complex).  Algorithm  A3  is  further  transformed  by  breaking  ys  into 
integer  and  fractional  parts:  ysj,  which  will  take  on  only  integer  values,  and  ysj-,  which  will  hold 
only  fractional  values.  Thus 

>s  =  ySi +  v  (2) 

0  <  ysf  <  1  (3) 

This  substitution  requires  that  the  incremental  step  0's  :  =  ys+[dy/dx J)  be  changed  to  add  the 
increment  to  the  fractional  part  (ysf)  and  then  test  whether  the  result  exceeds  1,  i.e.,  to  see  if  it  is 
no  longer  fractional. 

A4:  var  ysf:  exactreal:  dx,  dy,  xi,  ysi:  integer; 

ysi  :=  0;  ysf  :=  1/2; 
for  xi  :=  0  to  dx  do  begin 

(•  assert  ysi  +  ysf=yt+ 1/2  *) 
display(xi.ysi); 

if  ysf+fdy/dx]  >  1  then  begin 
ysi  :=  ysi  + 1; 
ysf  :  =  ysf  +  [dy/dx  -  1) 
end  else  begin 

ysf  :=  ysf + (dy/dx) 


end 
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4.  Substitution  of  variable  (simple).  Algorithm  A4  can  be  easily  transformed  into  the  Bresenham 
algorithm  by  replacing  the  use  of  ysj  with  that  of  a  variable  r: 

r  =  2  Jv  +  2(>'sy—  l)dx  (4) 

The  objectives  of  this  transformation  are  ( 1 )  to  change  the  comparison  in  the  inner  loop  to  a  sign 
check  (i.c.,  a  comparison  with  0).  and  (2)  to  eliminate  division  operations  by  scaling  by  2 dx. 
Making  the  appropriate  substitution  of  r  into  A4  yields  the  Bresenham  algorithm: 

A5:  var  dx.  dy.  xi.  ysi,  r.  integer; 

ysi  :=  0;  r  :=  2*dy  — dx; 
for  xi  :  =  0  to  dx  do  begin 

(•  assert  yi  vUt  »  ysf :  ysi  ♦  (lr  ♦  }dx  -  2dy)/2dx  *) 

display(xi.ysi); 
if  r  >  0  then  begin 
ysi  :  =  ysi  + 1 : 
r  :=  r-|2*dx-2*dyj 
end  else  begin 
r  :=  r  +  [2*dy] 
end 
end 

The  Bresenham  algorithm  is  ideal  for  implementation  in  hardware  or  microprocessors  with  limited 
arithmetic  power.  The  algorithm  requires  neither  division  nor  multiplication,  and  requires  no 
"floating-point"  approximations  because  all  variables  take  on  only  integer  values.  Moreover,  r  is 
not  required  to  hold  large  values.  Liquations  (3)  and  (4)  imply 

2dy-2dx  <  r<  2 dy  (5) 

If  0  <  dy  <  dx  <  2"— 1,  r  is  bounded  by 

-2"  +  1+2  <  r<  2,+1-I  (6) 

Thus  if  dx  and  dy  are  //-bit  positive  integers,  r  requires  at  most  n+2  bits  in  a  two’s  complement 
representation. 

Interpretation  of  r.  Ihc  value  of  r  is  Hated  to  the  vertical  error.  er  the  distance  from  the  pixel 
center  to  the  true  line.  Ihc  errors  will  be  identical  for  all  algorithms,  because  the  same  sequence  of 
points  is  generated.  When  display  is  called,  ev  =  )'sj~yr  Using  (1)  to  substitute  for  yf.  and  then 
(2)  to  substitute  for  we  have 

ev  =  L’jj- O'i-  1/2)  =  ysj—  (ysj+y’sj~  ^2) 

Cv  -  1/2  ~ysf 

Applying  transformation  (4)  yields 
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ev  =  -(r+dx—  2dy)/(2dx)  or 

r/(2dx)  =  -ev+(dy/dx- 1/2) 

The  value  r  is  thus  linearly  related  to  but  is  offset  by  1/2  due  to  the  loop's  initial  conditions, 
and  is  moreover  scaled  by  2 dx  to  require  only  integral  values  of  r. 

Summary:  Ah  of  the  algorithms  developed  in  this  section  compute  the  same  sequence  of  points 
(jr ,  y.)  that  approximate  the  true  line.  Mathematical  and  program  transformations  are  used  to 
derive  efficient  implementations. 


4.  The  digital  differential  analyzer  (I)DA) 

The  digital  differential  analy/cr  numerically  integrates  the  line  equation,  obtaining  x  =  {dx 
and  y  =  {dy.  'The  conventional  DDA  treats  both  coordinates  symmetrically.  The  numerical 
integration  requires  choosing  the  number  of  integration  steps,  as  shown  in  the  following  algorithm: 

A6:  var  x,y:  exactreal;  dx,  dy,  nsteps:  integer; 

x  :=  0;  y  :=  0; 

nsteps  :  =  some  number  >  dx  and  >  dy, 
for  i  :=  0  to  nsteps  do  begin 

display!  trunc(  x  +  [  1  /2]),trunc(y  +  [  1  /2])); 
x  :=  x  +  [dx/nstcps]; 
y  :  =  y +  [dy/nstcps) 
end 

This  algorithm  generates  exact  values  for  x  and  y  in  the  loop  because  "exact"  real  arithmetic  is 
used.  This  algorithm  may  not  produce  the  optimal  line,  in  the  sense  defined  in  Section  2,  because 
of  the  separate  rounding  in  x  and  y.  A  trivial  example  arises  if  nsteps  =  I :  only  the  pixels  at  the 
two  endpoints  of  the  line  will  be  displayed.  A  more  interesting  example  arises  when  dx  =  10,  dy 
~  8,  and  nsteps  is  chosen  to  be  40:  the  point  (2,  1)  is  displayed,  even  though  its  vertical  error  is 
-0.6. 

An  important  problem  with  the  DDA  is  the  choice  of  nsteps.  A  common  approach  is  to  choose 
nsteps  to  be  a  power  of  two  so  that  the  divisions  may  be  performed  simply  by  shifting  dx  and  dy: 
thus,  nsteps  =  2",  where  2"  >  dx.  Unfortunately,  this  causes  x  to  be  incremented  by  a  quantity 
less  than  unity,  so  that  more  than  one  pixel  in  a  column  may  be  illuminated.  Although  the  second 
pixel  in  a  column  can  be  omitted  by  a  suitable  test  in  the  loop,  it  is  harder  to  guarantee  that  the 
pixel  that  is  displayed  is  the  one  closest  to  the  line. 

Another  approach  is  the  "unit  increment"  DDA,  in  which  we  choose  nsteps  =  dx  so  that  x  will 
always  be  incremented  by  unity,  and  algorithm  A6  becomes  identical  to  A2.  which  generates  the 
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optimal  line.  By  this  route,  the  DDA  transforms  into  the  Bresenham  algorithm. 

It  is  important  to  note  that  a  common  approximation  to  the  unit  increment  DDA,  often  used  in 
hardware  implementations,  docs  not  generate  optimal  lines.  The  approximation  is  obtained  from  A4 
by  substituting  ysj  for  0-n)ysj-  in  order  to  introduce  an  integer  variable  yJ(.y  that  has  "sufficient" 
precision  to  represent  the  fractional  part  of  the  y  coordinate. 

A7:  var  dx.  dy,  xi.  ysi.  ysd:  integer; 

ysi  :=  0;  ysd  2n~l; 

for  xi  :=  0  to  dx  do  begin 
display(xi.ysi); 

if  ysd  +  [2n*(dy/dx)-  e)  >  2n  then  begin 
ysi  ysi  -t- 1; 

ysd  :=  ysd  +  [2n*(dy/dx)-e-2nj 
end  else  begin 

ysd  :=  ysd  +  (2"*(dy/dx)-e] 

end 

end 

The  above  algorithm  is  precise  only  if  e=0.  In  practice.  2 ”(dy/dx)-i  is  chosen  to  be  integral,  and 
e  is  the  error:  -1/2  <  e  <  1/2.  'Ihe  value  of  n  is  chosen  to  be  sufficiently  large  that  errors 
introduced  by  the  approximation  arc  acceptably  small.  If  n  is  chosen  to  be  the  smallest  integer  such 
that  2"  >  dx.  the  line  is  guaranteed  to  begin  and  end  at  the  proper  coordinates,  but  not  to  be 
optimal.  To  illustrate  the  non-optimality,  consider  dx  =  120,  dy  =  1,  n  =  7.  At  xi  =  62,  the 
algorithm  will  call  display (62,  0).  However,  the  point  (62,  1)  is  closer  to  the  true  line.  Ihere  is  a 
value  of  n,  2"  »  dx.  for  which  the  generated  line  will  be  optimal.  However,  this  value  may  be  quite 
high. 

The  hardware  implementation  of  the  unit-increment  DDA  is  simpler  than  algorithm  A7 
implies  (Figure  3).  On  each  iteration,  which  corresponds  to  a  clock  cycle  in  the  circuit,  a  new  value 
for  y is  computed  by  adding  2n(dy/dx)~  e,  >J(  will  be  incremented  if  the  sum  equals  or  exceeds 
2",  and  x,  will  always  increment.  The  same  idea  can  be  used  to  build  a  fixed-point  software 
implementation. 
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n  bits 
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Figure  3.  Hardware  implementation  of  the  unit-increment  DDA. 
Counters  compute  x(  and  y. .  An  adder  sums  yS(j  and  a  — 
2 n(dy/dx)-t.  On  each  cycle,  yS(j  is  loaded  with  the  new  sum,  x(-  is 
incremented,  and  •  is  incremented  if  the  high-order  bit  of  the 
adder  result  is  1. 


5.  An  n-step  algorithm 

Before  exploring  algorithms  that  exploit  parallelism,  we  shall  illustrate  the  transformation 

techniques  developed  in  the  previous  sections  by  deriving  an  algorithm  that  takes  horizontal  steps  of 
n  units  in  x.  Such  an  algorithm  will  generate  every  z/th  point  on  the  line.  We  start  with  an  obvious 
variant  of  Al: 

NT.  var  yt:  exactreal;  dx,  dy,  xi,  yi,  n:  integer; 

for  xi  ;=  0  do  dx  by  n  do  begin 

yt  ;  =  [dy/dx|*xi; 

yi  :=  trunc(yt+ 1/2); 

display(xi.yi) 
end 

Computing  yt  incrementally,  and  substituting  ys  =  yt+ 1/2,  we  have; 

N3:  var  ys:  exactreal:  dx,  dy,  xi,  yi,  n:  integer; 

ys  :=  1/2; 

for  xi  :  =  0  to  dx  by  a  do  begin 

yi  :=  trunc(ys); 

display(xi.yi); 
ys  :=  ys  +  [n*(dy/dx)) 
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When  ys  is  broken  into  integer  part  ysi  and  fractional  part  ysj-,  nUly/dx)  may  also  have  an  integer 
and  fractional  part.  Define  the  integer  part  s  so  that  0  <  n{dy/dx)-s  <  1;  the  "fractional"  part  is 
then  n(dy/dx)-s,  which  although  called  fractional,  may  actually  equal  1.  A  value  of  s  that  meets 
these  constraints  is  s  =  Ln(dy/dx)J.  Ihe  algorithm  becomes: 

N4:  var  ysf:  exactreal:  dx.  dy.  xi,  ysi.  n.  s:  integer; 

(•  assume  s  has  been  computed  *) 
ysi  :  =  0;  ysf  :=  1/2; 
for  xi  :  =  0  to  dx  by  n  do  begin 
display(xi,ysi); 

if  ysf-f  |n*(dy/dx)-s]  >  1  then  begin 
ysi  :=  ysi  +  (s+lj; 
ysf  :  =  ysf+(n*(dy/dx)-s— 1] 
end  else  begin 
ysi  :=  ysi  +  s; 
ysf  :=  ysf+|n*(dy/dx)-s] 
end 
end 

Ihe  next  step  is  to  apply  the  transformation  that  makes  a  "Brescnham-likc"  algorithm:  r  =  Indy  + 
2  (ysj~  1  -  s)dx. 


N5:  var  ysf:  exactreal:  dx,  dy,  xi,  ysi,  n.  s,  t:  integer; 

(•  assume  s  and  t  =  2nd)  -  2sdx  have  been  computed  *) 
r  :=  t— dx;  (*  ysf  =  1/2  implies  r  =  2ndy  + 2(1/2- 1 -s)dx  *) 

ysi  :=  0; 

for  xi  :=  0  to  dx  by  n  do  begin  "N5!oop” 

display(xi.ysi); 
if  r  >  0  then  begin 
ysi  :  =  ysi  +  [s+l]; 
r  :=  r-[2*dx-tj 

end  else  begin 
ysi  :=  ysi+s; 
r  :=  r+t 

end 

end  "N51oop" 

Note  that  this  algorithm  is  identical  to  A5  if  «=1.  s=0.  ihe  attentive  reader  will  question  what 
happens  if  dy-dx.  n=I.  Note  that  s  is  not  defined  to  be  Ln(dy/dx)d.  So  by  setting  s  =  0  in  this 
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case,  the  assumption  0  <  n(dy/dx)—s  <  1  is  not  violated.  1  he  other  possibility  for  dy=dx ,  j=  1, 
generates  the  same  points,  although  the  algorithm  is  then  not  identical  to  A5. 

A  minor  difficulty  with  N5  is  the  need  to  compute  s  and  i  -  2iuh  -  2sd\  Although  this  could 
be  done  with  multiply  and  divide  operations,  a  small  incremental  algorithm  can  he  used  to  compute 
s,  developed  using  the  same  principles  shown  in  A1-A5: 

var  nn:  exactrcal:  s.  i,  n:  integer; 

s  :  =  0.  nn  :  =  0; 

for  i  ;  -  0  to  n  - 1  do  begin 

if  rm  ■*- |d\/dx|  >  1  then  begin 

s  :  =  s+  1; 

rm  ;  =  rm  +  |dy/dx  - 1] 
end  else  rm  -  im-i-jdy/dx] 

end 

ITtis  program  is  transformed  by  substituting  t>  -  (nn  l),/t  +  </i  and  including  obvious  calculations 
for  I  into  the  following  prologue  to  algorithm  N5: 

N5p:  var  dx.  dy.  s.  t.  rp.  i.  n:  integer; 
begin  "NSprologuc" 
s  :  =  0;  t  :=  0; 
rp  :=  dy  — dx; 

for  i  :  =  0  to  n  •  1  do  begin 

;  (*  a*scrt  ildv/dvl  \-irp*d’  ri>|/dx  *1 

L  t  ;  =  t  +  dy; 

if  rp  >  t)  then  begin 

|  s  ;  =  s  a- 1 ; 

.  i  l  ;  -  t'-dx: 

rp  :  =■  rp  (dx  -  dy] 
end  else  rp  :  -  rp-i-dy 

end; 

t  ;=  t  +  t; 
end  ”N5pro(ogue" 

!  ft  is  important  to  remember  that  the  n  step  algorithm  generates  the  same  optimal  points  as  the 

j  Brcscnham  algorithm. 
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6.  Parallel  algorithms 

This  section  develops  line-drawing  algorithms  that  are  capable  of  generating  several  points  on  a 
line  in  parallel.  Ihese  algorithms  arc  useful  if  a  frame  buffer  can  update  several  pixels  in  one  cycle, 
or  if  lines  must  he  approximated  with  special  "characters"  |T  4].  The  transformations  illustrated  in 
the  preceding  sections  arc  extremely  useful  in  designing  these  algorithms.  They  allow  the  algorithm 
to  be  stated  in  conceptually  simple  terms  and  then  transformed  into  one  that  can  be  efficiently 
implemented  with  integer  arithmetic. 


6.1.  The  (n.n)  algorithm 

The  n-step  algorithm  developed  in  Section  5  is  the  basis  for  a  parallel  algorithm:  operate  n 
copies  of  the  procedure,  each  generating  points  spaced  n  units  apart:  hence  the  name  (n.n).  Hach 
copy  of  the  algorithm  is  phased  slightly  differently:  the  copy  with  phase—  0  generates  points  at  jk  =  0, 
n.  In.  .  .  .;  the  copy  with  phase=  I  generates  points  at  x=  1.  ;i+I,  2u+l.  .  .  .;  and  so  on.  Ibis 
technique  is  simply  expressed  as  (c.f.  A I  >: 

PI:  var  phase:  integer: 

for  phase  :  =  0  to  n  - 1  do  parbegin 

var  xi.yi:  integer:  yt:  exact  real.  (*  rhese  variables  are  duplicated  for  each  phase  *) 
for  xi  :=  0  + phase  to  dx  by  n  do  begin 
yt  :  =  [dy/dx]*xi; 
yi  :=  trunc(yt  +  (l/21); 
display(xi.yi) 
end 
pa  rend 

The  bracketing  parbegin  and  parend  mean  that  there  are  n  parallel  copies  of  the  inner  loop,  each 
operating  with  a  different  value  of  phase  and  with  separate  copies  of  the  local  variables  xi,  yi  and 
yt.  We  now  proceed  with  transformations  demonstrated  in  Sections  3-6.  lhe  inner  loop  is 
transformed  into  one  almost  identical  to  the  inner  loop  of  N5;  only  the  iteration  of  xi  is  different. 
The  initial  computation  for  ys  in  P2  requires  a  multiply/divide  which  is  transformed  into  a  loop 
executed  phase  times  to  compute  initial  values  for  yst  and  r.  This  loop  is  combined  with  the 
prologue  (N5p)  to  compute  values  for  s  and  t.  lhe  final  result  is  P2: 
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P2:  var  d*.  dy.  s.  t,  n.  rp.  i.  phase:  integer; 

begin  "N5prologue"  (*  Prologue  is  identical  lo  N5p  above  *) 
s  :=  0;  t  :=  0: 
rp  :=  dy— dx; 

for  i  :  =  0  to  n  - 1  do  begin 

(•  assert  Kdy/dx)  =  s-t(rp  +  dx-dy)/dx  •) 

t  :=  t+dy; 

if  rp  >  0  then  begin 

s  :=  s+1; 

t  ;=  t-dx; 
rp  ;=  rp-|dx-dy] 
end  else  rp  ;=  rp  +  dy 

end; 

t  :=  t+t; 
end  "N5prvlogue" 

for  phase  :=  0  to  n-1  do  parbegin 

var  xi.  ysi.  r.  i:  integer;  (*  This  variables  arc  duplicated  for  each  phase  *) 
r  :=  t-dx; 
ysi  ;  =  0; 
begin  "P2init" 

for  i  :=  0  to  phase -1  do 
if  r  >  |t-2*dy)  then  begin 
ysi  :  =  ysi  + 1 ; 
r  :=  r-[2*dx-2*dy] 
end  else  r  -  r+|2*dy] 
end  "P2init'’ 

for  xi  :=  0+ phase  to  dx  by  n  do  begin  "i  21oop" 

(•  Exactly  the  same  code  as  NSIoop.  above.  *) 
display(xi.ysi); 
if  r  >  0  then  begin 
ysi  :=  ysi  +  [s+1]; 
r  :=  r — [2*dx  —  t] 
end  else  begin 
ysi  :=  ysi+s; 
r  :=  r+t 
end 

end  "PIloop" 
pa rend 

Notice  that  the  loop  P2init  corresponding  to  the  initial  computations  for  ys  hears  a  strong 
resemblance  to  the  1  -step  Brcscnham  algorithm.  A5;  the  difference  arises  because  of  the  slightly 
different  expressions  for  r.  Thcic  is  no  particular  virtue  to  executing  the  loops  in  parallel— one  loop 
can  be  used  to  compute  initial  values  of  ysi  and  r  for  all  phases. 
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6.1  The  (l.n)  algorithm 

The  second  algorithm  capable  of  exploiting  parallelism  uses  the  //  step  algorithm  to  find  points 
on  the  line  at  //-unit  intervals  and  fills  points  in  between  with  a  "stroke."  The  n  pixels  in  each 
stroke  can  be  written  in  parallel.  This  technique  is  useful  when  lines  musi  he  approximated  with 
"characters"  because  a  raster  displa>  or  printer  is  controlled  b>  a  character  generator;  the  characters 
arc  simply  short  strokes. 

Ihc  algorithm  is  easily  derived  from  N5.  In  the  inner  loop,  the  test  on  r  determines  whether 
the  line  rises  by  s+1  or  s  units  for  a  move  of  /;  units  in  x.  If  the  line  rises  by  s  +  1  units,  a  stroke 
that  rises  s+1  units  in  n  is  drawn  from  the  current  (x  r)  point.  Ihc  stroke  is  determined  by  an 
index  /  that  gives  its  rise  in  v.  i  =  0.  1.  .  .  .  n.  The  strokes  mas  be  precomputed  using  the 
Bresenham  algorithm,  as  shown  in  Figure  4  for  /;- 8.  Note  that  each  stroke  has  only  n  points 
(jr  =  0,  1.  ...  n-  1).  but  that  the  rise  is  that  of  the  (n+  l)st  point  (x  =  n).  This  convention  is  adopted 
because  algorithm  N5  computes  the  rise  to  the  origin  of  the  next  stroke  rather  than  the  rise  to  the 
end  of  the  current  stroke. 
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Figure  4.  The  nine  different  strokes  for  /i~^.  The  left  column 
shows  rises  of  0  (bottom).  1,  2,  3,  and  4  (top).  The  right  column 
shows  rises  of  5  (bottom),  6,  7,  and  X  (top).  I  he  origin  of  a  stroke 
is  marked  with  a  +  and  the  origin  of  the  next  stroke  with  an  x. 
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In  order  to  draw  lines  of  arbitrary  length,  the  last  stroke  on  the  line  may  be  only  a  partial 
stroke.  The  standard  stroke  is  simply  truncated:  only  the  first  few  points  on  it  are  actually 
displayed.  Ihis  is  illustrated  by  the  procedure  DisplaySiroke.  which  accesses  an  array  .SVroArefex]  to 
find  the  y  coordinate  of  a  pixel  given  the  stroke  rise  t  and  the  x  coordinate  relative  to  the  beginning 
of  the  stroke. 

procedure  DisplayStrokc(originX,  originY,  rise.  maxX:  integer); 
var  x:  integer; 

for  x  :=  0  to  maxX  do  parbegin 

display(originX  +  x,  originY  +  Strokc{rise,x]) 

pa rend; 

Note  that  the  individual  pixels  of  the  stroke  arc  written  in  parallel. 

This  procedure  can  be  incorporated  into  N5  to  yield  the  complete  line-drawing  algorithm  Q. 
The  algorithm  is  shown  without  the  prologue  N5p: 

Q:  var  dx,  dy,  xi,  ysi.  s,  t.  r:  integer; 

(•  Insert  N5p  here  to  compute  s  and  t  -.  2ndy-2sdx  *) 
r  ;=  t-dx; 
ysi  :  =  0; 

for  xi  :=  0  to  dx  by  n  do  begin 
if  r  >  0  then  begin 

DisplayStrokcfxi.  ysi.  s+1.  min(n- l,dx-xi)) 
ysi  :=  ysi  +  (s+l]; 
r  :=  r-[2*dx  — t] 
end  else  begin 

DisplayStrokc(xi.  ysi.  s.  min(n- l.dx  —  xi)) 
ysi  :=  ysi  +  s; 
r  :=  r+t 
end 
end 

Algorithm  Q  has  several  advantages  over  P.  Ihc  setup  is  substantially  simpler,  as  are  the 
computations  performed  in  parallel.  The  scheme  is  very  similar  to  that  of  a  character  generator  in 
which  pre-computcd  patterns  arc  displayed:  the  strokes  play  the  role  of  characters.  It  differs  from 
many  character  generators  in  that  a  character  may  have  an  arbitrary  origin  on  the  screen  and  may 
be  partially  truncated. 

Ihe  chief  disadvantage  of  algorithm  Q  is  that  it  docs  not  generate  optimal  lines.  Although  the 
stroke  origins  lie  within  1/2  unit  of  the  true  line,  the  other  points  along  the  stroke  may  err  by  as 
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much  as  1  unit.  Iliis  property  arises  because  the  i  coordinate  of  a  pixel  is  the  sum  of  two 
independent  computations,  the  position  of  the  stroke  origin  and  the  position  of  the  pixel  withir  the 
stroke,  each  of  w  hich  may  make  an  error  of  i/2.  An  example  of  a  vertical  error  of  0.913  is  shown 
in  the  top  line  of  Figure  6.  at  x  -  18.  Another  way  to  see  the  non-optiinality  of  Q  is  to  observe 
that  although  only  a  single  stroke  is  displayed  for  each  distinct  rise  in  i.  there  are  actually  several 
different  strokes  with  the  same  rise  (Figure  5).  In  practice,  the  error  is  hardly  noticeable. 
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Figure  5.  Four  of  the  8  different  strokes  with  n  =  8  and  a  rise  of  1. 

Before  leaving  the  subject  of  stroke  selection,  we  should  mention  that  it  is  essential  to  have  the 
algorithm  P2  choose  from  two  strokes,  rather  than  merely  position  the  origin  of  a  single  stroke.  If  a 
single  stroke  is  used  for  an  entire  line,  the  maximum  deviation  from  the  optimal  line  may  be  greater 
than  1  or  the  line  may  have  gaps  or  non-monotonicities,  illustrated  in  Figure  6. 

F.vcn  though  algorithm  Q  produces  non-optimal  lines,  it  turns  out  that  the  endpoint  of  the  line 
is  always  exact.  Appendix  A  contains  a  proof  of  this  fact.  Fxhaustive  simulations  of  algorithm  Q 
for  all  lines  of  length  1024  or  less  have  verified  that  endpoints  are  always  computed  correctly. 

7.  Conclusion 

This  paper  began  by  showing  how  simple  mathematical  and  program  transformations  could  be 
used  to  transform  an  obvious  line-drawing  method  based  on  analytic  geometry  into  an  efficient  and 
exact  algorithm  that  requires  only  integer  arithmetic.  Iliese  methods  help  persuade  us  that  the 
algorithm  is  correct  without  recourse  to  complex  geometric  constructions  such  as  those  in  [1],  Ihc 
techniques  arc  examples  of  routine  program  transformations  that  should  be  a  commonplace  activity 
in  program  design  and  implementation. 
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Figure  6.  l  ines  illustrating  gaps  and  non-monotonicitics.  The  top 
line  (dx  =  23.  dy= 18)  is  drawn  with  three  strokes  with  «  =  8,  which 
leave  a  gap.  The  small  dots  show  the  optimal  line.  Ihe  bottom 
line  (dx-  23,  dy=  5)  shows  a  non-monotonicity. 


The  main  reason  for  applying  these  techniques  is  to  extend  line-drawing  algorithms  to  exploit 
parallel  activities.  Although  only  two  parallel  schemes  arc  explored  in  Section  6.  one  can  imagine 
many  more.  Ihe  difficulty  of  developing  such  algorithms  is  subtantially  reduced  by  using  the 
program  transformations. 
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Ihis  paper  grew  from  attempts  to  write  very  fast  line-drawing  microcode  for  the  "8x8  display." 
designed  by  Ivan  Sutherland  and  the  author.  Satish  Gupta  devoted  considerable  coding  effort  to 
this  display  and  to  simulations  of  the  (lot)  method.  The  proof  in  Appendix  A  is  due  to  Mike 
SprciUcr  of  Caltech. 
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Appendix  A 

We  demonstrate  that  the  (l.n)  algorithm  terminates  at  the  proper  endpoint  Idx.dy)  Assume  dx 
>  dv  >  0.  let  2  =  ax  +  by  be  the  measure  of  distance  from  the  line,  and  further  let  a  -  2 dy.  b  = 
-  2dx.  Define  v(  to  be  the  y  coordinate  o.  the  pixel  displayed  at  x  =  i:  z,  is  the  distance  from  this 
pixel  to  the  true  line. 

The  Bresenham  algorithm  that  generates  the  origin  for  a  stroke  will  guarantee  that  |zj  < 
-b/2.  The  points  .*  =  ni  +  j  for  /  =  1  U)  n  arc  members  of  one  of  the  two  strokes  that  represent 
a  line  of  slope  s/n.  where  s  is  the  vertical  distance  from  the  origin  of  the  stroke  to  the  origin  of  the 
next  stroke  (i.e..  s  =  >  -  ynl).  If  these  pixels  arc  generated  by  a  Bresenham  algorithm  aiming 

at  a  line  of  slope  s/n.  then  we  will  have  ~  z„,)l  <  —b/2.  for  0  <  j  <  n—  1. 

We  consider  two  cases:  first,  that  the  expression  is  positive,  and  second,  that  it  is  negative. 

1.  We  have  zn  <  (//n)(zm+n  -  z„)-b/2.  From  the  triangle  inequality,  we  also  have 
|z  _  t  |  <  -b  However,  the  equality  ease  never  occurs — if  it  did.  the  slope  of  the  line 

' ni+n  nr  — 

would  be  an  integer  multiple  of  \/n  and  the  rm  would  be  zero  for  all  /.  So  we  now  have:  | zm>(|  - 
Zm I  <  -b.  For  1  <  ;  <  n/2.  this  yields  zm+j  <  -  b. 

2.  The  negative  ease,  by  similar  argument,  gives  znj+ ,  >  b  for  1  <  j  <  n/2. 

Both  cases  together  give  |zm+  |  <  -b  for  1  <  ./  <  n/2.  By  a  similar  argument  approaching 
from  the  other  side  (i.e..  x  =  ni  +  n.  ni  +  n  —  1,  .  .  .).  we  obtain  |-m_y|  <  ~  b  for  1  <  j  <  n/2. 
Both  forward  and  backward  approaches  together  give  \zmJrJ\  <  -b  for  1  <  <  n. 

When  x  =  dx.  at  the  endpoint  of  line,  we  must  have  <  —  b.  which,  together  with  the  fact 
that  z  must  be  a  multiple  of  b  forces  zJx  =  0.  Ihcrefore.  the  last  point  lies  exactly  on  the  line. 


